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Abstract 

We explore numerically the evolution of a collapsing spherical shell 
of charged, massless scalar field. We obtain an external Reissner- 
Nordstrom space-time, and an inner space-time that is bounded by 
a singularity on the Cauchy Horizon. We compare these results with 
previous analysis and discuss some of the numerical problems encoun- 
tered. 

1 Introduction 

The most general analytic, stationary, and asymptotically flat solution of 
the Einstein equations is the Kerr-Newman metric, which is characterized 
by three parameters - mass, charge and angular momentum. The no-hair 
theorem [Tj suggests that at asymptotically late times after the collapse of a 
massive body to a black hole, the external metric of the black hole relaxes 
to the Kerr-Newman metric. However, the final state of the inner region is 
not generally determined and has been the subject of research for some time. 
Also, it is of interest, especially to gravitational wave astronomers, to predict 
the specific details of the relaxation of the "hair" . 
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Figure 1 : A Penrose diagram for the maximally extended Reissner-Nordstrom 
space-time. The Cauchy Horizon is marked as r = r_, The Event Horizon 
as r = r + . Taken from [2]. 

This appearance of a time-like singularity and a "chain" of asymptotically 
flat regions in the maximally extended Reissner-Nordstrom metric (Fig.(jT|)) 
(and also in the Kerr metric) opened up the discussion of the existence and 
stability of "wormholes". The Reissner-Nordstrom metric is problematic as 
a realistic black hole picture because the part of space-time that lies in the 
future of the Cauchy Horizon is unpredictable by a set of Cauchy initial 
conditions given on a space-like surface in the asymptotically flat region. 
This happens because events in the region beyond the Cauchy Horizon are 
in causal connection to the time-like singularity. 
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Consequently we would expect that a in a dynamic collapse the space-time 
will not extend beyond the Cauchy Horizon, but would rather be bounded 
there by a singularity of some sort. This scenario was predicted by Penrose 
[Sj as early as 1968. The physical mechanism leading to this singularity is 
known as the "Mass Inflation" scenario 0], and is connected to the infinite 
blueshift of in-falling radiation on the Cauchy Horizon. 

Previous results [3 EH IZj concerning the stability of the cauchy horizon 
were acquired either by linear perturbation theory or by numerically evolving 
a neutral perturbation collapsing on a preexisting Reissner-Nordstrom black 
hole [SI EI] • The latter model is relatively easy to solve numerically since 
it does not contain charges and dynamical electromagnetic fields, but on the 
other hand it does not give the complete picture of the evolution of the black 
hole and the Cauchy Horizon beginning with a flat space time. Previous 
work dealing with a fully dynamical charged collapse model EH ITi] 

focused mainly on the behavior near the origin and outside the horizon, i.e 
at relatively large values of retarded time but moderate values of advanced 
time, on the order 10M. These works indicated that in a charged collapse 
the Cauchy Horizon is replace by a weak, null singular segment connected 
to a space-like singularity at the origin (Fig.(j2j)). In this work we follow the 
evolution of a charged matter shell that collapses to form a black hole, and 
determine the final structure of the inner metric, and the Cauchy Horizon at 
large values of advanced time, deep inside the Black Hole interior. 

In section 2, we formulate our basic equations, starting from the gauge 
invariant Lagrangian of the complex scalar field and casting these equations 
into the double-null coordinate system. In section 3 we discuss the numerical 
scheme, it's accuracy and the difficulties that arise in this kind of calculation. 
In section 4 we discuss the numerical results and their implications regarding 
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Figure 2: A Penrose diagram for a charged collapse scenario. The dashed 
lines indicate singularities along the origin and Cauchy Horizon. Taken from 

® 

several aspects of black hole physics, such as Mass Inflation, Cosmic Cen- 
sorship, Wave tails and the No-Hair theorem, and the nature of the Cauchy 
Horizon singularity. 

2 Formulation of the problem 

2.1 The equations 

We begin with the gauge-invariant Lagrangian for the complex scalar field if) 
and the electromagnetic gauge field jT3] : 

— 77?^ — 1 

% = -(tp. a + teA^)g ab (^ b - ieA^) - - —F ab F cd g ac g bd , (1) 

where we define = A u -^ — A mv . We are dealing with a massless field so 
we will omit the middle term from now on. From this lagrangian we can 
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derive the Euler-Lagrange equations for the scalar and electromagnetic field, 
and also the Energy-Momentum tensor to be used in the Einstein equations. 
The wave equations for the scalar field are: 

i{j ]ab g ab + ieA a g ab {2^, b + ieA^) + ieA a , b g ab ^ = 0, (2) 

and its complex conjugate. The Maxwell equations take the form: 

7^ F ab; C g bc + ieip{^a - ieA a ^) - ie^a + ieA a ij) = 0. (3) 

The Energy-Momentum tensor given by: 

Tab = -(V>;a^ ;b + Vv0 ; &) + ^{-^- a ieA b ^ + ^. b ieA a ^ + ^. a ieA b ^ - ^. b ieA a ^) + 
^F ac F bd g cd + e 2 A a A b ^ + ±S? 9ab . (4) 

We choose the (2+2) spherically symmetric, double null coordinate sys- 
tem. The line element in these coordinates can be written as: 

ds 2 — —a(u, v) 2 dudv + r(u, v) 2 dVt 2 , (5) 

with r being the area coordinate. Note that any coordinate gauge of the form 
u — > f(u),v — > g(v) will preserve the null character of u and v. We will fix 
this gauge freedom later when discussing the initial conditions. 

We need also to fix the gauge freedom of the electromagnetic field A^. In 
our symmetry and coordinate choice the only non- vanishing components of F 
are F uv = —F vu . This means that only A u and A v may be non-zero. We can 
eliminate one of these using the electromagnetic field gauge freedom, A^ — > 

+ 3> ;/1 , where <E> is an arbitrary scalar function. We choose to eliminate A v 
by taking $ = — J A v dv, and so we are left with A^ = (A u , 0, 0, 0). We shall 
call this quantity a = A u from now on. 
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We now proceed to write our equations in an explicit form. We begin with 
the Maxwell equations. Since we reduced the potential to one component 
using the gauge freedom, we need only one equation for this field. Eq.(|3j) 
is a vector equation with two non-vanishing components u and v, and we 
may choose either one to evolve the electromagnetic field. We choose the 
v component, and will use the u component later when checking charge 
conservation. The v component then becomes (we are using from now on 
the convention Z a = J^r): 

T 2 Cb 1 

( — ir-)v + -ier 2 K{^ v - ipt{j v ) = 0. (6) 

We define: 

2r 2 n 

Q = — 5-, (7) 

cr 

and thus separate the 2 nd order equation for a into two simpler 1 st order 
equations: 

a 2 q 

a " = w 

q v = 2mear 2 (ipifj v — ij>if) v )- (8) 

The function q(u, v) defined above is the amount of charge within (i.e. at 
smaller radii than) a sphere of radius r(u,v) on some space-like hypersurface 
that contains (u,v). 

The mutually independent elements of G^ v read: 

. 2oj 1 ,7 > 1/ (XT*. 



G U u 2 ■ 



G vv 

G U v 



l u' u L *' uu 

ar 

^2(y v T v CXT VV 

ar 

a 2 + 4r u r v + 4rr v 
2r~ 2 



ee = 4 , (yj 

cr 



Combining Eq.fjUJ) and the Einstein tensor G^ v we arrive at the field equa- 
tions: 

a v 

r vv -2r v h Anr%% = 

a 

r uu ~ 2r„— + \-Kr{i) u ^ u + iea(ipi) u - $ij) u ) + eV^) = 
a 

a 2 a 2 q 2 
rr uv + r u r v + — - — = 

— V 1 + — + —tt + ^(tfuifiv + $ v i/j u ) + 2Tiiea{^ v - i)^ v ) = (10) 

a or r 4r 4 

Finally we evaluate the wave equation and arrive at the form: 

a 2 q 

ri>uv + r u ib v + r v ip u + iearip v + iear v ip + ie——ijj = 0. (11) 

4r 

We are now in a position to write down the full set of equations to be 
solved. For this purpose we introduce some new notations: 



47T0 

^ — $ u , z = s v 

f = r u ] g = r v (12) 
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Rearranging and substituting we arrive at the following set of equations: 

a 2 a 2 a 
El: rr uv + f 9 + ---^ = 

P , x f9 a 2 a 2 q 2 1 1. _ 

El : {lna) uv - — + -t-t- + -{wz + wz) + -iea(sz - sz) = 

r z 4r z 2r 4 2 2 

Oi 

CI : o„ — 2— <? + rzz = 
a 
cx 

C2 : f u — 2—f + riuw + ieaisw — sw) + e 2 a 2 ss 
a 

2 

Ml: a„-|| = 

M2 : ^ — -ier 2 (sz — sz) = 

a 2 q 

S: rs uv + fz + gw + iearz + ieags + ie— — s = 0. (13) 

4r 



2.2 The initial conditions 

Having specified the equations to be solved we turn to discuss the formulation 
of the initial conditions. The physical situation we wish to describe is the 
gravitational collapse of a shell of in-falling charged matter. We choose the 
domain of integration to be a rectangle in the u-v plane (see section 3.2, 
where we describe the numerical scheme). We can expect the Event Horizon 
to be inside the domain of integration, but not the Cauchy Horizon since 
it is located at infinite null coordinate v (We are not using Kruskal-like 
coordinates), we can However approach the Cauchy Horizon asymptotically 
at large values of v. 

When we come to specify the metric functions on the initial hypersurface, 
we must fix the coordinate gauge freedom mentioned in section 2.1. In a 
Minkowsky space-time we usually choose u = t — r,v = t + r. Space-time 
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Figure 3: Illustration of the domain of integration in relation to the expected 
structure of the space-time. The Cauchy Horizon cannot be actually included 
in the domain because it is in null infinity. Taken from jTHj 

around a gravitating spherical shell is flat in two regions - inside the shell 
and at asymptotically large radii. Since we have one constraint equation, 
CI, that relates r,a and the scalar field on the initial hypersurface, we can 
choose two of these functions arbitrarily. We would like to take r to be 
linear with v in order to reflect the fact that we start our integration far 
from the Event Horizon , where the metric is nearly flat. The fact that this 
is only an approximation will be pronounced by the deviation of a from its 
flat-space value. We can further specify the shape of the matter shell on the 
initial hypersurface, and integrate CI to obtain a on that surface. Using the 



9 



conventional flat-space dependence of r on u and v, 

n = — g— , (14) 

we conclude that = 0.5, and given an initial distribution of the scalar field 
that we will soon elaborate on, we get from CI the following ODE for a: 

— = rzz. (15) 
a 

We choose for the initial field distribution a compactly supported function, 
which is basically one half period of a cosine. Studying equation M2 we see 
that in order to have a non-vanishing charge distribution we must use for the 
real and imaginary part of the field two pulses offset in v by some amount. 
The exact form we choose is: 

Vi > v > 1>2 



Re[si 



Ar( v 2 ' Y vi < v < v 2 

v[ > v > v' 2 



Ai C z^' ) 2 v[<v<v' 2 . 
Having specified the two fields r and s on the initial hypersurface, we can 
analytically derive z and g (which are simply g — \ and z = ||), and nu- 
merically integrate the remaining quantities /, a, a, q and w using equations 
El, CI, Ml, M2 and S respectively. 

We also need to specify the boundary conditions on the ray v — v {. Since 
the metric on the boundary, which is always inside the collapsing shell, is 
flat, we can select the conventional flat-space values for the fields there: 

Vi - Ui 

r = 

2 

a = 1 

s = z = q = A = 
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(17) 



In selecting these values for the metric functions on the boundary we specify 
the physical meaning of the metric functions: The choice of a implies that 
^Y^- is the proper time for an observer located at the origin, while r is the 
proper surface area of a sphere of a radius r(u, v). 

For monitoring the mass content in our space-time we use the mass func- 
tion, 

2 

nr* Ci T 1 T* 

m(u,v) = -(l+ q - + 4^f), (18) 
which represents the total mass in a sphere of radius r(u, v). 

3 The numerical scheme 

3.1 Domain of integration 

The main consideration in choosing the domain of integration is whether 
or not to include the origin. The origin will be important when discussing 
critical phenomena ^7], however we are only interested in effects that occur 
at large values of v and at a finite radius, so we can choose the domain to 
be a rectangle as shown in Fig.fjHJ). The advantage of this approach is that 
the origin, while being regular, is a coordinate singularity, and this would 
force us to make a series expansion of the fields near the origin in order to 
maintain numerical precision. 

3.2 Numerical Algorithm 

We construct a 2-dimensional grid in the v,u space and integrate in increasing 
v and u direction. Since our equal-coordinate surfaces lie on the light-cone, 
there are no Courant-like stability limits on the step size. We implement a 
Runge-Kutta like scheme by making a trial step to the middle of the cell 
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and using this information to make a 2 nd order full step. We have found this 
scheme to be simple and efficient. 

We use equation El to evolve r,f and g. E2 is used to evolve a, and S 
evolves s,z and w. Finally Ml and M2 evolve a and q respectively. Note, that 
this scheme is a Free Evolution scheme, since only the dynamical equations 
are used for evolving the space-time. The constraint equations are used for 
monitoring the accuracy. 

3.3 Adaptive grid refinement 

A major problem in black hole numerics is the behavior near the Event Hori- 
zon . Although the double-null coordinates we have chosen ensure regularity 
of the fields even as we cross the horizon, there is nevertheless a fundamen- 
tal difficulty in following the evolution numerically because of the physical 
behavior near the Horizon. Consider two outgoing null rays starting at the 
origin, slightly above and below the Event Horizon. One is destined to es- 
cape to infinity while the other will remain trapped inside the Horizon. This 
means that regardless of how close the two rays were initialii, their distance 
will diverge as their advanced time v grows. The numerical implication is 
that the metric function f — r u diverges along the Event Horizon, and r 
becomes "discontinuous" asymptotically. If we want to maintain a fixed rel- 
ative change in the area coordinate r, we have to make the grid denser as 
we approach the horizon. The step size has to decrease with the distance 
form the horizon. This decay turns out to be exponential, and the exponent 
increases as we increase v max . It would be best to change both Au and Av 
as we approach the horizon and reach large values of v, respectively, but 
we chose in this work to change only Au because while this is significantly 
simpler to implement, it gives reasonable results. The choice of having the 
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Figure 4: Adaptive grid refinement. The left figure features the development 
of u along a typical grid, showing how most of the computational effort is 
concentrated at the neighborhood of the Horizon. The figure on the left is 
the difference Au, showing that Am changes by many orders of magnitude. 

refinement algorithm maintain — constant, instead of some other local indi- 
cator of accuracy, is not trivial and was made after trial and error indicated 
that it was the best strategy, and also because r is always finite in the re- 
gion of interest so it provides a good scale for measuring error. This scheme 
dictates, for relatively high values of v max ~ 100M, a condensation of about 
10 orders of magnitude in Au in the vicinity of the Event Horizon relative 
to the initial value of Au. This is illustrative of the difficulty in numerically 
crossing the Event Horizon at large values of v. Fig.(|4*|) shows an example of 
the evolution of Au on a typical space-time. 
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3.4 Numerical tests 



Our first test would have been a comparison of the numerical results with 
a known analytical solution. Unfortunately, there are no suitable analytical 
solutions. Therefore, we will have to check our code by other means. We 
begin with a test of convergence, i.e. we verify that the various fields converge 
to some value in the expected order of convergence. We evolve the same 
initial conditions on three grids, g%, gi and g% which have fixed grid densities 
(in both u and v directions) d\, 2d\ and Ad\. We need three grids because we 
do not know the exact solution, so we can compare only the relative error. In 
Fig.(j2J) we plot, for example, the functions Re[ipi — ^2] and 4Re[ip2 ~ ^3]- We 
expect, and verify, that the two curves be nearly but not exactly the same, 
because higher order terms also have a small contribution to the error. The 
data in this figure is from a single row(w = const) but it is representative of 
the entire grid. 

Although convergence is crucial, it still does not guarantee that the results 
are correct, because if there is an error in the equations, either in the orig- 
inal equations or as they are implemented in the code, the scheme will still 
converge but to the wrong solution. We are in a unique situation where the 
constraint equations provide us with a measure of consistency. The Einstein 
equations preserve the constraints, in the sense that if the initial conditions 
satisfy the constraint equations, then they will be satisfied also under evo- 
lution of the dynamical equations. This is true analytically, but since our 
solution is numerical we cannot solve exactly either the constraint equations 
on the initial hypersurface or the dynamical equations through the evolution. 
Therefore we can only demand that the error in the constraint will converge 
to zero as we make our solution more accurate by condensing the grid, i.e. 

stepsize— >0 

C > . This is verified in Fig.fjHJ). 
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Figure 5: Re[ipi — ip 2 ] & 4:Re[ifj 2 — ^3] for the real part of the scalar field on 
one of the outgoing rays. The two curves are almost indistinguishable, which 
shows that the scheme converges nicely to 2 n d order 



Next we check that our code reproduces known features, such as the 
location of the horizons, mass and charge conservation etc. Comparing the 
location of the horizons to the Reissner-Nordstrom solution's r± = M ± 
a/ M 2 — q 2 we find a correspondence of about one part in 10 5 between the 
computed and expected radii. 

Finally we observe that global mass and charge are conserved. Actually 
this is not exactly the case since some of the field is scattered by the gravita- 
tional and electromagnetic potential as the shell approaches its gravitational 
radius, but this is in most cases (when Q < M) an insignificant amount and 
the conservation laws seem to hold. Fig. (f7jl shows the mass M and the charge 
Q at the last grid-point vs. u. Also shown is Q' which is Q compensated for 
the charge lost at the last grid point due to outgoing flux: 

Q'(u) = Q(u)+ r 2 (u',v max )Ju(u',v max )du', (19) 

15 
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Figure 6: The maximum constraint error along outgoing rays max v (Ci), 
max v (4:C2) and max v {lQC^) as a function of u on the same 3 grids as in 
Fig- - The three curves are again almost indistinguishable, showing that 
the constraint violation displays 2 nd order convergence to zero. The initial 
error in the constraint grows with u because of unavoidable accumulation of 
numerical errors. 

which follows from the conservation equation for the electromagnetic current: 

J* = 0. (20) 

It should be noted that the mass inside the domain of integration is not 
generally conserved, since the reflected waves carry off some mass, as well 
as charge. The subject of mass conservation is more delicate than charge 
conservation since the total mass contains also the Energy- Momentum of the 
gravitational field, which doesn't enter into the Noether current related to 
the Energy-Momentum tensor. We will refrain from dealing with this topic 
here. 
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Figure 7: Global charge and mass conservation. Upper figure: Q is the charge 
q(u, v), Q' is the charge, compensated for the outgoing current density at the 
last grid point (outgoing flux). Lower figure: the relative change 1 — ^ and 
1 — which is four orders of magnitude smaller. 

3.5 Numerical error 

To summarize, we can obtain a measure of the relative numerical errors 
from the tests we discussed. The first test is to use the convergence test 
and compare the difference in the fields between different grids. There is a 
difficulty here because this test can only be performed without using adaptive 
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Figure 8: Right panel: The metric function r vs. v for several outgoing null 
rays, on 3 different grids #1,2,3 with grid spacing 4h, 2h and h respectively. 
The three grids are compared to a solution extrapolated from them, taking 
into account the 2 nd order convergence. This solution is much more accurate 
than the three others and can be taken to as the "accurate" solution for 
the purpose of comparison. The 4 curves are almost indistinguishable. Left 
panel: detail of the leftmost ray 

grids. This substantially limits the precision that can be attained in this 
test, but as Fig.® illustrates, even with this limitation we get good results: 
~ 10~ 3 maximum relative error. As we indicate in the next paragraph, the 
error will be lower by an order of magnitude or more if we employ the grid 
refinement algorithm. 

The next test is charge conservation. We can see in Fig.® that the 
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charge conservation error is less than 10~ 4 for the most dense grid. We can 
also see that the error converges to zero. These results are taken from refined 
grids, because in this case we are not comparing rows against each other so 
we don't need to match rows exactly in u value. This test was done with 
v max — 20M. In this case the same calculation on an unrefined grid gave an 
error larger by an order of magnitude. This ratio will be larger as we increase 
v max , making the use of grid refinement more and more crucial. 
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Figure 9: 1 — It- for three different grids vs. u. The dashed line indicates the 
position of the EH. The legend indicates the relative step size used in each 
grid. 

The last test we perform is constraint violation. We can split CI in 
Eq.(|12jl into two terms and compare their relative error: 

A = g v -2% (21) 
a 

B = -rzz. (22) 

I M I _g| 

Now we can take ^ + # to be the relative error. However there is a problem 
since \A\ + \B\ has no definite scale and it often intersects zero, so the error 
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Figure 10: The relative error in the constraint, max l 
The three curves are for the same 3 grids as in Fig.(jHJ). 
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behaves very badly. Therefore we will compare (A — B) in each u = const 
segment to the r.m.s. of (A + B) along the same segment. The results are 
given in Fig. fjlOj) . and show a maximum error level of ~ 10~ 5 in the finest 
grid used. 

4 Results 



4.1 Formation of the black hole 

We begin Fig- lfTTj) with a null-ray diagram giving an overview of the space 
time. The diagram shows the formation of an Event Horizon and then a 
Cauchy Horizon at asymptotically large v. Looking closely at the null rays 
approaching the Cauchy Horizon we see that it is not stable, as the asymp- 
totic (v — > oo) value of r decreases as u increases. Finally at large enough 
values of u the rays reach the origin (r = 0), which signals the appearance 
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Figure 11: Null rays. The top figure shows the entire domain of integration, 
and the three classes of rays: type-1 rays reach null-infinity, type-2 approach 
the Cauchy Horizon and type-3 terminate at the singular origin. The dashed 
line signifies the approximate location of the Apparent Horizon, the two 
lower figures focus on the vicinities of rays approaching the Cauchy Horizon 
(left) and the Event Horizon (right). 

of a space-like singularity. 
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4.2 charge Dynamics 

The Cosmic Censorship conjecture [TJ] states that all singularities in nature 
are contained within an Event Horizon , so there are no "Naked singularities" . 
It follows from this that it is impossible to overcharge a black hole (i.e. 
to increase its charge-to-mass ratio above 1) since there is no black hole 
solution with this ratio. It is interesting to check if this holds in our collapse 
scenario by increasing the initial charge of the collapsing shell. It can be 
seen in Fig.(|12j) that, for a certain initial configuration, as we increase the 
coupling constant e, and with it the initial charge to mass ratio, the final 
black hole charge to mass ratio reaches a maximal value smaller than 1 
and then decreases. The physical mechanism that causes this behavior is 
the electrostatic repulsion of the outer parts of the shell, which increases 
relative to the gravitational pull they experience as the charge increases, and 
causes them to be reflected back towards large radii at the early stages of the 




Figure 12: The charge to mass ration on the Event Horizon for different 
initial conditions. 
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collapse and radiate part of the initial charge away. This can be understood 
as a scattering process, where at each instant in time the shell encounters 
some potential barrier with certain reflection and transmission coefficients. 
As long as the charge-to-mass ratio will be greater than unity, no horizon will 
form, and charged matter will be able to escape. Ordinarily the charge and 
mass become equal some finite time before the shell crosses it's gravitational 
radius and a Horizon forms, so the final charge-to-mass ratio decreases further 
and ends up smaller than unity. 

This process can be seen by inspecting the movement of charges on the 
computational grid. The outgoing current density, in our coordinates, is 
proportional to the u component of the 4-current J M . Note that this is a 
matter of interpretation since normally the charge density would be the time- 
like component of the 4-current, and the current density would be the space- 
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like component. In our coordinate system u and v are both null, but we treat 

v and u as the "spatial" and "temporal" coordinates because we are dealing 

with an in-falling shell. Fig. (|T3~j) shows a contour of this quantity for the case 
Qo = i 7 

Mo ' ■ 

4.3 Wave tails 

Next we observe the mechanism behind the No-Hair theorem, mentioned 
in chapter 1. We inquire how the black hole radiates away the "hair", i.e. 
any feature of the collapsing matter except it's charge, mass and angular 
momentum. We want to compare this decay to theoretical results predicting 
a certain power law decay of the tails at late times. Note that we perform 
this analysis along the Event Horizon, which is a null hypersurface with a 
constant radius. Hod&Piran ^3] predicted in this case the following behavior 
on the Horizon: 

^ ~ e if ? v~ p ; (3 = 1 + v/(2/ + l) 2 -4(eQ) 2 , (23) 

where I is the multipole moment of the scalar perturbation. In our case 
I = because we deal with spherical symmetry. The first term in Eq. (j23|) 
is an oscillation term with angular frequency — . The second term is a 
power-law decay of the field magnitude. Typically the field magnitude on the 
horizon begins to decay with Quasi-Normal ringing (exponentially decaying 
oscillations) which decays exponentially, Then the power-law tail sets in and 
continues to asymptotically large advanced time v where it dominates the 
behavior of the field. This can be clearly seen in Fig.(|14jl. 

The above expression for the oscillation and power law index were checked 
over a range of the parameter eQ, by changing e and leaving all other pa- 
rameters constant. The results in Fig. (jl5j) show good correspondence with 
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Figure 14: A typical configuration of the scalar field amplitude on the 
horizon, showing Quasi-Normal ringing followed by a power-law decay (the 
"tail"). The thick line is the amplitude of the field, the thin line is the 
absolute value of the real part of the field. 
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both terms in Eq. (J23|) . 




Figure 15: Upper figure: The real part of the power law index for spherical 
charged perturbations versus the theoretical prediction. Lower figure: The 
phase oscillation frequency versus the theoretical value, 

4.4 The Cauchy Horizon 

We now turn to examine the behavior of the space-time near the Cauchy 
Horizon. First we try to verify that null rays of type-2 actually reach an 
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Figure 16: The behaviour of r v along a null ray approaching the Cauchy 
Horizon . The power-law decay can be seen clearly at late times. 
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Figure 17: The local power-law index of r v at the approach to the Cauchy 
Horizon . r v behaves like v~ 2 at large v, signifying that r itself approaches a 
constant value as The noise at large v is only a numerical artifact. 
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asymptotically fixed radius, rather than assuming a course that would ulti- 
mately bring them, at very large values of v, to the origin. In order to do 
this we must extend our domain of integration, so that we can observe the 
decay of r v over at least one order of magnitude and establish the asymptotic 
behavior. In Fig. (jl7j) and Fig. f|17|) such an analysis is made, showing that 
r v ~ v~ 2 . This confirms the existence of a Cauchy Horizon. 

The next thing we want to know is wether the Cauchy Horizon is singu- 
lar, and if so how much. One indication of a singularity is an exponential 
divergence of the mass function. This can be seen in Fig. lfTHJ) . for a typical 
space-time. This divergence is dominated by an exponential decay of a with 
v, which makes the mass function diverge exponentially, because it depends 
on a as a' 2 . This can be interpreted physically as an infinite blue-shift of 
in-falling radiation at the Cauchy Horizon, since a is also gu if we shift to 
a time-space (1+3) coordinate frame. The theoretical prediction for mass- 
inflation is m ~ e KV , where k is the surface gravity on the inner horizon. For 
a static Reissner- Nordstrom black hole this is given by [H] Eq.(J2H), and it 
agrees numerically to about 10 — 20% with the measured exponent. 



The singularity also manifests itself in the Ricci curvature scalar, R. The 
Ricci scalar is a purely geometrical entity, which can be written only with the 
metric functions and their derivatives. However, we can use our dynamical 
equations to express this quantity through the physical fields in the problem. 
We arrive at this form: 



k = 




(24) 




a 



(25) 
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Figure 18: Mass inflation near the CH. Each curve is for a constant u, u is 
increasing from lower to higher lines. 

R can be seen to diverge exponentially on the horizon while oscillating at 
some definite frequency (see Fig. p^jl ). This behavior indicates the existence 
of some kind of curvature singularity, but does not tell us if it is strong or 
weak, i.e. if the tidal forces felt by an observer crossing it are infinite or 
finite. The divergence in R is caused by the exponential decay of a, and the 
oscillation is caused by the oscillatory behavior of the scalar field. We can 
observe the frequency of this oscillation by looking at Ra 2 which contains 
the oscillating term without the exponential divergence. Fig.(|20p shows this 
oscillation and it's Fourier transform. A strong oscillation with period 0.5M 
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Figure 19: The Ricci scalar on an outgoing null ray intersecting the Cauchy 
Horizon. The cusps indicate points where R changes sign. 

is accompanied by a weaker component with period 0.6M, and an oscillation 
with "beat" is produced. 

It is interesting to note here, that as the curvature radius at the Horizon 
of a black hole is on the order of magnitude of it's gravitational radius, which 
is ~ 10 4 m for a solar mass black hole, and the planck scale is ~ 10~ 35 m, the 
Ricci curvature can grow by approximately 40 orders of magnitude before it 
reaches Planckian scales. This means that our analysis brings us well into 
the Plackian regime, where classical GR begins to fail. 

5 Summary and conclusions 

We have developed a numerical code that solves the coupled Einstein-Maxwell 
equations in a dynamical collapse situation. This has enabled us to explore 
phenomena involved in black hole formation that was previously handled by 



30 





Figure 20: The oscillatory part of the Ricci scalar. Top panel shows an os- 
cillation with a weak "beat" . The bottom panel shows the Fourier transform 
of the curve in the top panel. 

analytic or numerical perturbative analysis. The difficulty in this kind of 
numerical code is in maintaining small numerical errors in spite of extreme 
gradients near the Event Horizon. We have solved this problem by using a 
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non-uniform grid that covers the difficult areas with dense grid points. The 
scheme was shown to converge and give results accurate to better than than 
1 part in 10 3 , even in areas deep inside the black hole. 

We have shown that in a dynamical collapse of charged matter, some of 
the charge is radiated away because of electrostatic repulsion and scattering 
on the gravitational potential. This is in accord with the Cosmic Censorship 
conjecture which forbids the charge in the Black Hole to surpass it's mass. 
We also observed the radiative tails that are left on the Horizon after the 
black hole is formed. These were shown to have the well known structure 
of an initial decaying oscillation (Quasi-Normal ringing) followed by a power 
law decay that continues asymptotically. 

Finally we ventured deep inside the black hole interior and examined the 
properties of the inner Horizon. We have found that before being completely 
destroyed and turning into a strong space-like (i.e. Schwarzschild-like) sin- 
gularity it behaves as a weak, null singularity. Since the Cauchy Horizon is 
located at v = oo in our scheme we cannot reach it by numerical evolution. 
However physically it is reached in a finite proper time by an in-falling ob- 
server (This is true since a 2 (u,v)dt which is the proper time differential for 
an observer at (u, v) decays exponentially as v — > oo, giving a finite lapse 
of proper time until reaching the Cauchy Horizon). The weakness of the 
singularity thus leaves open the question of the traversability of the "Kerr 
Tunnel" , making it unclear wether it is physically possible for matter to cross 
the CH into a another asymptotically flat region. 
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